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The diversity of neuron models used in contemporary theoretical neuroscience to investigate 
specific properties of covariances in the spiking activity raises the question how these models 
relate to each other. In particular it is hard to distinguish between generic properties of 
covariances and peculiarities due to the abstracted model. Here we present a unified view 
on pairwise covariances in recurrent networks in the irregular regime. We consider the binary 
neuron model, the leaky integrate-and-fire model, and the Hawkes process. We show that 
linear approximation maps each of these models to either of two classes of linear rate models, 
including the Ornstein-Uhlenbeck process as a special case. The distinction between both 
classes is the location of additive noise in the rate dynamics, which is located on the output 
side for spiking models and on the input side for the binary model. Both classes allow closed 
form solutions for the covariance. For output noise it separates into an echo term and a term 
^1 due to correlated input. The unified framework enables us to transfer results between models. 

qv For example, we generalize the binary model and the Hawkes process to the situation with 

£■ — synaptic conduction delays and simplify derivations for established results. Our approach 

I* is applicable to general network structures and suitable for the calculation of population 

averages. The derived averages are exact for fixed out-degree network architectures and 
approximate for fixed in-degree. We demonstrate how taking into account fluctuations in 
the linearization procedure increases the accuracy of the effective theory and we explain 
the class dependent differences between covariances in the time and the frequency domain. 
Finally we show that the oscillatory instability emerging in networks of integrate-and-fire 
models with delayed inhibitory feedback is a model-invariant feature: the same structure of 
poles in the complex frequency plane determines the population power spectra. 
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Correlations, linear response, Hawkes process, leaky integrate-and-fire model, binary neuron, 
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1. Introduction 

The meaning of correlated neural activity for the processing and representation of informa- 
tion in cortical networks is still not understood, but evidence for a pivotal role of correlations 
increases (recently reviewed in Cohen & Kohn, 2011). Different studies have shown that 
correlations can either decrease (Zohary et al., 1994) or increase (Sompolinsky ct al., 2001) 
the signal to noise ratio of population signals, depending on the readout mechanism. The ar- 
chitecture of cortical networks is dominated by convergent and divergent connections among 
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the neurons (Braitenberg & Schiiz, 1991) causing correlated neuronal activity by common 
input from shared afferent neurons in addition to direct connections between pairs of neurons 
and common external signals. It has been shown that correlated activity can faithfully prop- 
agate through convergent-divergent feed forward structures, such as synfire chains (Abeles, 
1991; Diesmann ct al., 1999), a potential mechanism to convey signals in the brain. Cor- 
related firing was also proposed as a key to the solution of the binding problem (von der 
Malsburg, 1981; Bicncnstock, 1995; Singer, 1999), an idea that has been discussed controver- 
sially (Shadlcn & Movshon, 1999). Independent of a direct functional role of correlations in 
cortical processing, the covariance function between the spiking activity of a pair of neurons 
contains the information about time intervals between spikes. Changes of synaptic coupling, 
mediated by spike-timing dependent synaptic plasticity (STDP, Markram et al., 1997; Bi 
& Poo, 1999), are hence sensitive to correlations. Understanding covariances in spiking 
networks is thus a prerequisite to investigate the evolution of synapses in plastic networks 
(Burkitt ct al., 2007; Gilson et al., 2009, 2010). 

On the other side, there is ubiquitous experimental evidence of correlated spike events 
in biological neural networks, going back to early reports on multi-unit recordings in cat 
auditory cortex (Perkcl et al., 1967; Gerstcin & Perkcl, 1969), the observation of closely 
time-locked spikes appearing at behaviorally relevant points in time (Maldonado et al., 
2004; Kilavik et al., 2009) and collective oscillations in cortex (recently reviewed in Buzsaki 
& Wang, 2012). 

The existing theories explaining correlated activity use a multitude of different neuron 
models. Hawkes (1971) developed the theory of covariances for linear spiking Poisson neurons 
(Hawkes processes). Ginzburg & Sompolinsky (1994) presented the approach of linearization 
to treat fluctuations around the point of stationary activity and to obtain the covariances for 
networks of non-linear binary neurons. This formal concept of linearization allowed Brunei 
& Hakim (1999) and Brunei (2000) to explain fast collective gamma oscillations in networks 
of spiking leaky integrate-and-fire (LIF) neurons. Correlations in feed-forward networks of 
leaky integrate-and-fire models are studied in Morcno-Bote & Parga (2006) , exact analytical 
solutions for such network architectures are given in Roscnbaum & Josic (2011) for the case 
of stochastic random walk models, and threshold crossing neuron models are considered in 
Tchumatchcnko et al. (2010) and Burak ct al. (2009). Covariances in structured networks are 
investigated for Hawkes processes (Pernice et al., 2011), and in linear approximation for LIF 
(Pernice et al., 2012) and exponential integrate-and-fire neurons (Trousdale et al., 2012). 
The latter three works employ an expansion of the propagator (time evolution operator) 
in terms of the order of interaction. Finally Buice et al. (2009) investigate higher order 
cumulants of the joint activity in networks of binary model neurons. 

The diversity of the employed models brings up the question which features of covariances 
are generic properties of recurrent networks and which are particular to a certain model. 
Currently it is unclear how different neuron models relate to each other and whether and 
how results obtained with one model carry over to another. In this work we present a 
unified theoretical view on pairwise correlations in recurrent networks in the asynchronous 
and collective oscillatory regime, approximating the response of different models to linear 
order. 

The remainder of this article is organized as follows. In the first part of our results in 
"Covariance structure of noisy rate models" we investigate the activity and the structure of 
covariance functions for two versions of linear rate models (LRM); one with input the other 
with output noise. If the activity relaxes exponentially after application of a short perturba- 
tion, both models coincide with the Ornstein-Uhlenbeck process (OUP). We mainly consider 
the latter case, although most results hold for arbitrary kernel functions. We extend the 
analytical solutions for the covariances in networks of OUP (Risken, 1996) to the neuroscien- 
tifically important case of synaptic conduction delays. Solutions are derived first for general 
forms of connectivity in "Solution of the convolution equation with input noise" for input 
noise and in "Solution of convolution equation with output noise" for output noise. After an- 
alyzing the spectral properties of the dynamics in the frequency domain in "Spectrum of the 
dynamics" , identifying poles of the propagators and their relation to collective oscillations 
in neuronal networks, we show in "Population-averaged covariances" how to obtain pairwise 
averaged covariances in homogeneous Erdos-Renyi random networks. We explain in detail 
the use of the residue theorem to perform the Fourier back-transformation of covariance 
functions to the time domain in "Fourier back transformation" for general connectivity and 
in "Explicit expression for the population averaged cross covariance in the time domain" for 
averaged covariance functions in random networks, which allows us to obtain explicit results 
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Figure 1: Mapping different descriptions of neuronal dynamics to linear rate models (LRM). 
The arrows indicate analytical methods which enable a mapping from the original spiking 
(leaky integrate-and-fire model, Hawkes model) or binary neuron dynamics to the analyt- 
ically more tractable linear rate models. Depending on the original dynamics (spiking or 
binary) the resulting LRM contains an additive noise component x either on the output side 
(left) or on the input side (right). 

and to discuss class dependent features of covariance functions. 

In the second part of our results in "Binary neurons" , "Hawkes processes" , 
and "Leaky integrate-and-fire neurons" we consider the mapping of different neuronal dy- 
namics on either of the two flavors of the linear rate models discussed in the first part. 
The mapping procedure is qualitatively the same for all dynamics as illustrated in Figure 1: 
Starting from the dynamic equations of the respective model, we first determine the working 
point described in terms of the mean activity in the network. For unstructured homoge- 
neous random networks this amounts to a mean-field description in terms of the population 
averaged activity (i.e. firing rate in spiking models). In the next step, a linearization of the 
dynamical equations is performed around this working point. We explain how fluctuations 
can be considered in the linearization procedure to improve its accuracy and we show how 
the effective linear dynamics maps to the LRM. We illustrate the results throughout by 
a quantitative comparison of the analytical results to direct numerical simulations of the 
original non-linear dynamics. The appendices "Implementation of noisy rate models" , 
"Implementation of binary neurons in a spiking simulator code" , and 

"Implementation of Hawkes neurons in a spiking simulator code" describe the model imple- 
mentations and are modules of our long-term collaborative project to provide the technology 
for neural systems simulations (Gewaltig & Diesmann, 2007). 

2. Covariance structure of noisy rate models 

2.1. Definition of models. Let us consider a network of linear model neurons, each 
characterized by a continuous fluctuating rate r and connections from neuron j to neuron 
i given by the element u>ij of the connectivity matrix w. We assume that the response of 
neuron i to input can be described by a linear kernel h so that the activity in the network 
fulfills 

r(t) = ft(o) * [wr(o - d) + bx(o)](t), (I) 

where f(o — d) denotes the function / shifted by the delay d, x is an uncorrelated noise with 



(Xi(t)) 



0. 



(xi(s)xj(t)) = 8 i:j 8(s - t)p 2 , 



(2) 



e. g. a Gaussian white noise and (/ * g)(t) = /_ f(t — t')g(t') dt' is the convolution. With 
the particular choice b = w8(o — d)* we obtain 

r(t) = [h(o) * w(r(o - d) + x(o - d))](t). (3) 
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We call the dynamics (3) the linear noisy rate model (LRM) with noise applied to output, 
as the sum r + x appears on the right hand side. Alternatively, choosing b = 1 we define 
the model with input noise as 

r(t) = h{o) * [wr(o -d) + x(o)](t). (4) 

Hence, equations (3) and (4) are special cases of (1). In the following we consider the 
particular case of an exponential kernel 

h(s) = -6(s)e~ s / T , (5) 

where 8 denotes the Heaviside function, 6(t) — 1 for t > 0, else. Applying to (1) the 
operator O — t-^ + 1 which has h as a Green's function (i.e. Oh — S) we get 

t— r(t)+r(t)=wr(t-d)+hx(t), (6) 

at 

which is the equation describing a set of delay coupled Ornstein-Uhlenbeck-processes (OUP) 
with input or output noise for b = 1 or b = w5(o — d)*, respectively. We use this represen- 
tation in "Binary neurons" to show the correspondence to networks of binary neurons. 

2.2. Solution of the convolution equation with input noise. The solution for the 
system with input noise obtained from the definition (4) after Fourier transformation is 

R = ff d wR + HX, (7) 

where the delay is consumed by the kernel function hd(s) = -8(s — d)e~( s ~ d ^/ T . We use 
capital letters throughout the text to denote objects in the Fourier domain and lower case 
letters for objects in the time domain. Solved for R = (1 — H^w)^ 1 HX the covariancc 
function of r in the Fourier domain is found with the Wiener-Khinchin theorem (Gardiner, 
2004) as (R(w)R T (— w)), also called the cross spectrum 

C(w) = (R(cj)R t (-o;)) (8) 

= (1 - H d (u;)vr)- l H(uj)(X{uj)X T {-uj))H(-uj)(l - Hdi-Lu)^)- 1 
= (ff d (w)- 1 -w)- 1 D( J ff d (- w )- 1 -w T )- 1 , 

where we introduced the matrix D = (X(w)X (— ui)}. From the second to the third line we 
used the fact that the non-delayed kernels H(uj) can be replaced by delayed kernels H^uj) 
and that the corresponding phase factors e lujd and e~ lud cancel each other. If x is a vector 
of pairwise uncorrelated noise, D is a diagonal matrix and needs to be chosen accordingly in 
order for the cross spectrum (8) to coincide (up to non-linear effects) with the cross spectrum 
of a network of binary neurons, as described in 
"Equivalence of binary neurons and Ornstein-Uhlcnbcck processes" . 

2.3. Solution of convolution equation with output noise. For the system with 
output noise we consider the quantity y$ = r, + Xi as the dynamic variable representing the 
activity of neuron i and aim to determine pairwise correlations. It is easy to get from (3) 
after Fourier transformation 

R = ff d w(R + X), (9) 

which can be solved for R = (1 — iJ^w^^iJ^wX in order to determine the Fourier transform 
of Y as 

Y = R + X = (l-J? d w)- 1 X. (10) 

The cross spectrum hence follows as 

C(«) - (Y(c)Y T (-c)) (11) 

= (1 - J ff d ( W )w)- 1 (X( W )X T (-^))(l - H d (-uj)w T )- 1 
- (1 - J ff d ( W )w)- 1 D(l - H d (-io)^ T )-\ 
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with D = (X(w)X T (— u))). D is a diagonal matrix with the i-th diagonal entry p\. For the 
correspondence to spiking models D must be chosen appropriately, as discussed in "Hawkes 
processes" and "Leaky integrate-and-fire neurons" for Hawkes processes and leaky integrate- 
and-fire neurons, respectively. 

2.4. Spectrum of the dynamics. For both linear rate dynamics, with output and 
with input noise, the cross spectrum C(w) has poles at certain frequencies w in the complex 
plane. These poles are defined by the zeros of &et(Hd{i>j)~ l — w) and the corresponding term 
with the opposite sign of ui. The zeros of det(iJrf(w) _1 — w) are solutions of the equation 



H d {u)~ l = (1 + iujT)e iud = L 



where Lj is the j-th eigenvalue of w. The same set of poles arises from (1) when solving for 
R. For d > and the exponential kernel (5), the poles can be expressed as 

z k {L j ) = ---W k {L ] -e^), (12) 

t a t 

where Wk is the fe-th of the infinitely many branches of the Lambert-W function (Corless 
et al., 1996). For vanishing synaptic delay d = there is obviously only one solution for 
every Lj given by z = =?(Lj - 1). 

Given the same parameters d, w, r, the pole structures of the cross spectra of both 
systems (8) and (11) are identical, since the former can be obtained from the latter by 
multiplication with (H^^H^—u)) -1 = {H{u)H{— w)) _1 , which has no poles. The only 
exception causing a different pole structure for the two models is the existence of an eigen- 
value Lj = of the connectivity matrix w, corresponding to a pole z(0) = -. However, 
this pole corresponds to an exponential decay of the covariance for input noise in the time 
domain and hence does not contribute to oscillations. For output noise, the multiplication 
with the term (H(u)H(-uj))^ 1 , vanishing at u> = -, cancels this pole in the covariance. 
Consequently both dynamics exhibit similar oscillations. A typical spectrum of poles for a 
negative eigenvalue Lj < is shown in Figure 2B,D. 

2.5. Population-averaged covariances. Often it is desirable to consider not the whole 
covariance matrix but averages over subpopulations of pairs of neurons. For instance the 
average over the whole network would result in a single scalar value. Separately averaging 
pairs, distinguishing excitatory and inhibitory neuron populations, yields a 2 by 2 matrix of 
covariances. For these simpler objects closed form solutions can be obtained, which already 
contain some useful information and show important features of the network. Averaged 
covariances are also useful for comparison with simulations and experimental results. 

In the following we consider a recurrent random network of N e = N excitatory and 
Ni = 7-/V inhibitory neurons with synaptic weight w for excitatory and —gw for inhibitory 
synapses. The probability p determines the existence of a connection between two randomly 
chosen neurons. We study the dynamics averaged over the two subpopulations by intro- 
ducing the quantities r a = jj-J2je a r j an( i noise terms x a = j$-'52j£a x j f° r a ^ {£>Z}; 
indices X and £ stand for inhibitory and excitatory neurons and corresponding quantities. 
Calculating the average local input -/V^ 1 Xwea w jk r k to a neuron of type a, we obtain 

c'EE^^ ( 13 ) 

j £ a k 

= N^ipNaiv yVfc- pN a gw 2J r k ) 

fcef kex 

= pwN(r £ -7flr z ), 

where, from the second to the third line we used the fact that in expectation a given neuron 
k has pN a targets into the population a. The reduction to the averaged system in (13) is 
exact if in every column k in Wj k there are exactly K non-zero elements for j € £ and 
^K for j £ X, which is the case for networks with fixed out-degree (number of outgoing 
connections of a neuron to the neurons of a particular type is kept constant), as noted 
earlier (TetzlafT et al., 2012). For fixed in-degree (number of connections to a neuron coming 
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Figure 2: Pole structure determines dynamics. Autocovariance of the population activ- 
ity (A,C) measured in p 2 /r and its Fourier transform called power spectrum (B,D) of the 
rate models with output noise (dots) (3) and input noise (diagonal crosses) (4) for delays 
d = 3 ms (A,B) an d d — 1 ms (C,D) delay. Black symbols show averages over the exci- 
tatory population activity and gray symbols over the inhibitory activity obtained by direct 
simulation. Light gray curves show theoretical predictions for the spectrum (20) and the 
covariance (21) for output noise and the spectrum (17) and the covariance (18) for input 
noise. Black crosses (12) in B, D denote the locations of the poles of the cross spectra - 
with the real parts corresponding to the damping (vertical axis), and the imaginary parts 
to oscillation frequencies (horizontal axis). The detailed parameters for this and following 
figures are given in "Parameters of simulations" . 



in from the neurons of a particular type is kept constant) the substitution of 



jEa 



with 



r a is an additional approximation, which could be considered as an average over possible 
realizations of the random connectivity. In both cases the effective population-averaged 
connectivity matrix M turns out to be 



M = Kw 



-19 
'19 



(14) 



with K — pN. So the averaged activities fulfill the same equations (3) and (4) with the non- 
averaged quantities r, x, and w replaced by their averaged counterparts r = (r£,rz) T , x = 
(x£,xx) T , and M. The population averaged activities r a are directly related to the block- 

C££ CSX 



wise averaged covariance matrix c = 
With 



cis cxi 



withc o6 = JV- 1 Jvr 1 E< e aE 



£>„, =N7 1 n: 




jEb Ci i 



(15) 



=8 ab N a /N 2 aP 2 = 8 ab N- l p 2 

N- 1 

(7JV)- 
(8) and their general solutions also hold for the block-wise averaged covariance matrices. 

The covariance matrices separately averaged over pairs of excitatory, inhibitory or mixed 
pairs are shown in Figure 2 for both linear rate dynamics (3) and (4). (Parameters for 



we replace D by D = p 2 



_i I and c by c so that the same equations (11) and 
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Figure 3: Limits of the theory for fixed in-degree and fixed out-degree. Autocovariance 
(A) and covariance (B) in random networks with fixed in-degrcc (dots) and fixed out-degree 
(crosses). Simulation results for C££, C£x, and c%x are shown in dark gray, black and light 
gray, respectively for synaptic weight w = 0.011 far from bifurcation. For larger synaptic 
weight io = 0.018 close to bifurcation (see text), egg is also shown in D for fixed in-degree 
(dark gray dots) and for fixed out-degree (black dots) . Corresponding theoretical predictions 
for the autocovariance (34) (A) and the covariance (18) (B,D) are plotted as light gray curves 
throughout. The set of eigenvalues is shown as black dots in panel C for the smaller weight. 
The gray circle denotes the spectral radius wy/ Np(l — p)(l + ^g 2 ) (Rajan & Abbott, 2006; 
Kriener et al., 2008) confining the set of eigenvalues for the larger weight. The small filled 
gray circle and the triangle show the effective eigenvalues L of the averaged systems for small 
and large weight, respectively. 



all simulations presented in this article are collected in "Parameters of simulations", the 
implementation of linear rate models is described in "Implementation of noisy rate models"). 
The poles of both models shown in Figure 2B are given by (12) and coincide with the peaks 
in the cross spectra (8) and (11) for output and input noise, respectively. The results of 
direct simulation and the theoretical prediction are shown for two different delays, with the 
longer delay leading to stronger oscillations. 

Figure 3C shows the distribution of eigenvalues in the complex plane for two random con- 
nectivity matrices with different synaptic amplitudes w. The model exhibits a bifurcation, 
if at least one eigenvalue assumes a positive real part. For fixed out-degree the averaging 
procedure (13) is exact, reflected by the precise agreement of theory and simulation in Fig- 
ure 3D. For fixed in-degree, the averaging procedure (13) is an approximation, which is good 
only for parameters far from the bifurcation. Even in this regime still small deviations of the 
theory from the simulation results are visible in Figure 3B. In the proximity of a bifurcation, 
the appearance of long living modes causes large fluctuations. These weakly damped modes 
appearing in one particular realization of the connectivity matrix are not represented after 
the replacement of the full matrix w by the average M over matrix realizations. The eigen- 
value spectrum of the connectivity matrix provides an alternative way to understand the 
deviations. By the averaging the set of N eigenvalues of the connectivity matrix is replaced 
with two eigenvalues of the reduced matrix M, one of which is zero due to identical rows 
of M. The eigenvalue spectrum of the full matrix is illustrated in Figure 3C. Even if the 
eigenvalue(s) L M of M are far from the bifurcation region (corresponding to 3(z(L M )) > 0) 
some eigenvalues L w of the full connectivity matrix in the vicinity of the bifurcation region 
may still have an imaginary part becoming negative and the system can feel their influence, 
shown in Figure 3D. 
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2.6. Fourier back transformation. Although the cross spectral matrices (8) and (11) 
for both dynamics look similar in the Fourier domain, the procedures for back transformation 
differ in detail. In both cases, the Fourier integral along the real w-axis can be extended to 
a closed integration contour by a semi-circle with infinite radius centered at in the appro- 
priately chosen half-plane. The half-plane needs to be selected such that the contribution of 
the integration along the semi-circle vanishes. By employing the residue theorem (Bronstcin 
et al., 1999) the integral can be replaced by a sum over residua of the poles encircled by 
the contour. For a general covariance matrix we only need to calculate c(i) for t > 0, as for 
t < the solution can be found by symmetry c(— t) — c T (— t). 

For input noise it is possible to close the contour in the upper half-plane where the 
integrand C(w) e lut vanishes for \u\ — > oo for all t > 0, as |C^(w)| decays as M~ 2 . This can 
be seen from (8), because the highest order of H^ oc oo appearing in det(H^ — w) is equal 
to the dimensionality N of w (N — 2 for M), and in det(adjugate matrix ijoi H^ 1 — w) it 
is N - 1 (i = j) or N - 2 (i ^ j). So {(Hj 1 - w) _1 | is proportional to |w|~ 1 |e~* a " i | and 
|C(w)| oc |w|" 2 for large |w|. 

For the case of output noise (11) C(w) can be obtained from the C(ui) for input noise 
(8) multiplied with {Hd{oj)Hd{—ijj))~ l ~ |cli| 2 for large \ui\. The multiplication with this 
factor changes the asymptotic behavior of the integrand, which therefore contains terms 
converging to a constant value and terms decaying like |w| _1 for |w| — > oo. These terms 
result in non-vanishing integrals over the semicircle in the upper half-plane and have to be 
considered separately. To this end we rewrite (11) as 

C(w) =((1 - ff d (w)w)- 1 H d (w)w + l)T)(w T H d (-u J )(l - Hdi-u)^)- 1 + 1) (16) 

= (1 - H d {u)w)- l H d (w)viTiw T H d (-w){l - H d (-u)w T )- 1 

+{l-H d (uj)w)- 1 H d {cj)yvT> 
+Dw T fl d (- W )(l - H d {-uj)^ T )- 1 
+D, 

and find the constant term D which turns into a (5-function in the time domain. The first 
term in the second line of (16) decays like \uj\~ 2 and can be transformed just as C(w) for 
input noise closing the contour in the upper half-plane. The second and third term are the 
transposed complex conjugates of each other, because of the dependence of H on — u> instead 
of w, and require a special consideration. Multiplied by e under the Fourier integral, the 
first term is proportional to H d e lut ~ UJ - 1 e lu >( t - d ) an d vanishes faster than |w| _1 for large \u\ 
in the upper half-plane for t > d and in the lower half plane for t < d. For the second term 
the half planes are interchanged. The application of the residue theorem requires closing 
the integration contour in the half-plane where the integral over the semi-circle vanishes 
faster than |cj| _1 . For w = M and in the general case of a stable dynamics all poles of 
the first term are in the upper half-plane Q(zk(Lj)) > 0, and have no contribution to c(i) 
for t < d. For the second term the same is true for t > — d; these terms correspond to the 
jumps of c(i) after one delay, caused by the effect of the sending neuron arriving at the other 
neurons in the system after one synaptic delay. These terms correspond to the response of 
the system to the impulse of the sending neuron - hence we call them "echo terms" in the 
following (Helias ct al., 2013). The presence of such discontinuous jumps at time points 
d and —d in the case of output noise is reflected in the convolution of hw with D in the 
time domain in (37). For input noise the absence of discontinuities can be inferred from the 
absence of such terms in (33), where the derivative of the correlation function is equal to 
the sum of finite terms. The first summand in (16) corresponds to the covariance evoked by 
fluctuations propagating through the system originating from the same neuron and we call 
it "correlated input term" . In the system with input noise a similar separation into effective 
echo and correlated input terms can be performed. We obtain the correlated input term 
as the covariance in an auxiliary population without outgoing connections and echo terms 
as the difference between the full covariance between neurons within the network and the 
correlated input term. The echo and correlated input terms for models corresponding to 
linear rate models with either input or output noise are illustrated in Figure 7. 

2.7. Explicit expression for the population averaged cross covariance in the 
time domain. We obtain the population averaged cross spectrum in a recurrent random 
network of Ornstcin-Uhlenbeck processes with input noise by inserting the averaged con- 
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nectivity matrix w = M (14) into (8). The explicit expression for the covariance function 
follows by taking into account all (both) eigenvalues of M with values and L = Kw(l — jg). 
The detailed derivation of the results presented in this section are documented in "Explicit 
calculation of the expression for the population averaged cross covariance in time domain" . 
The expression for the cross spectrum (8) takes the form 

C(u>) = mf{-u>)(l + Kw(™ -J°^H d (u)\l)(l + Kw( 2 9 ig \)H d {-u 

(17) 

where we introduced /(w) = (Hdiuj)^ 1 — L)^ 1 as a short hand. Sorting the terms by 
their dependence on w, introducing the functions $i(cj), . . . , $4(0;) for this dependence, and 
(pi(t),. . . , ip4,(t) for the corresponding functions in the time domain, the covariance in the 
time domain c(t) = ^ J C(w)e l "*do; takes the form 

c(*) = D^i(t) 

+ *v(7 -_T)°(-; a )«<«>■ 

The previous expression is valid for arbitrary D. In simulations presented in this article 
we consider identical marginal input statistics for all neurons. In this case the averaged 
activities for excitatory and inhibitory neurons are the same, so we can insert the special 
form of D given in (15), which results in 

c(f, =£(i 7°- )*"<" (18) 

+ tH™ -;->H + £*»(™ -7-H 

The time-dependent functions (pi, . . . , (pi are the same in both cases. Using the residue 
theorem Vi (t) = ^ f+™ $;(w)e ia,i dw = iE zep oics of a>, Res($i, z) e"* for i ^ they can 
be expressed as a sum over the poles Zk(L) given by (12) and the pole z = - of Hd(u>)- 
At w = %(£) the residue of /(w) is Res(/,w = Zfe(L)) = (idL + ire lujd ) , the residue of 
H e i(uj) at z = - is —-e d l T 1 so that the explicit forms of <^i, ... ,<pt follow as 

Vl (t)= J2 *Res(/, W )/(- W )e^ 

u=z h (L) 

^(t)= X! iRes(/, W )/(-a;)F d ( W ) e - t + ^^/(-)/(-i) 

* — ' r r r 

u=z fc (i) 

p 3 (*) = E «Res(/,u;)/(-u;)F (J (-u;)e^ (19) 

uj=z h (L) 

Mt) = E ^es(/, W )/(- W ) J ff d (c)i/ d (- W ) e - t + £^/(i)/(-I). 
The corresponding expression for C(w) for output noise is obtained by multiplying (17) with 



C(«) - H^\u)Hj\-w)f(w)f{-u) (20) 

x (1 + fto ( Y "I 7 / ) ffd(u;))D(l + Kw ( 2! jg \ ) H d {-w)) 
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which, after Fourier transform, provides the expression for c(t) in the time domain for t ^ 
c(t) = MDM T ^(() + MD^ (t) + D6(t) 

-kv£(i+™">(! ;)„w + a-4(; :j)«,m+£(; £)«*>. 

(21) 

As in (18), the first line holds for arbitrary D, and the second for D given by (15), valid if 
the firing rates are homogeneous, ipi is defined as before, and 

(po(t) = 0(t - d) J2 {dL + re^y 1 e luJt (22) 

u>=z k (L) 

vanishes for t < d. All matrix elements of the first term in (21) are identical. Therefore 
all elements of c(i) are equal for < |t| < d. Both rows of the matrix in front of ipo are 
identical, so for t > the off diagonal term c%e coincides with C£ £ and C£x with c%x an( A Yice 
versa for t < 0. The solution (18) is visualized in Figure 5 and the solution (21) in Figure 6. 

3. Binary neurons 

In the following sections we study, in turn, the binary neuron model, the Hawkes model 
and the leaky integrate-and-fire model and show how they can be mapped to one of the two 
OUPs; either the one with input or the one with output noise, so that the explicit solutions 
(18) and (21) for the covariances derived in the previous section can be applied. In the 
present section, we start with the binary neuron model (Ginzburg & Sompolinsky, 1994; 
Buicc ct al., 2009). 

Following Ginzburg & Sompolinsky (1994) the state of the network of N binary model 
neurons is described by a binary vector n 6 {0, 1}^ and each neuron is updated at inde- 
pendently drawn time points with exponentially distributed intervals of mean duration r. 
This stochastic update constitutes a source of noise in the system. Given the i-th neuron is 
updated, the probability to end in the up-state (n, = 1) is determined by the gain function 
Fi(n) which depends on the activity n of all other neurons. The probability to end in the 
down state (n, = 0) is 1 — .Fj(n). Here we implemented the binary model in the NEST 
simulator (Gewaltig & Diesmann, 2007) as described in "Implementation of binary neurons 
in a spiking simulator code" . Such systems have been considered earlier (Ginzburg & Som- 
polinsky, 1994; Buicc ct al., 2009), and here we follow the notation employed in the latter 
work. In the following we collect results that have been derived in these works and refer the 
reader to these publications for the details of the derivations. The zero-time lag covariance 
function is defined as Cij(t) — (rii(i)nj(i)} — ai(t)a,j(t), with the expectation value () taken 
over different realizations of the stochastic dynamics. Here a(£) = (ai(t), . . . , ajy(t)) T is the 
vector of mean activities a,i(t) — (rii(t)). Cij(t) fulfills the differential equation 



d 



r^Ciiit) = -2a s (t) + (K(t) - a*(t)) Ji(n)) + ((n^t) - a^F^n)). 



In the stationary state, the correlation therefore fulfills 

di = i<(n i -a,)F i (n)) + i((r k -a i )F,(n)}. (23) 

The time lagged covariance cy(t, s) = (ni(t)nj(s)) — ai(t)a,j(s) fulfills for t > s the differential 
equation 



d 

~dt 



T—Cij(t,s) = -c iJ (t,s) + ((n j (s)-a j (s))F l (n 1 t)). (24) 



This equation is also true for i = j, the auto covariance. The term ((fij(s) — aj(s))Fi(n,t)} 
has a simple interpretation. It measures the influence of a fluctuation of neuron j at time s 
around its mean value on the gain of neuron i at time t (Ginzburg & Sompolinsky, 1994). 
We now assume a particular form for the coupling between neurons 



N 



Fi(n,t) = ^{3 i Tn{t-d))=(j>C£ d J ik nk{t-d)), (25) 



fe=i 
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Figure 4: Alternative linearizations of the binary neuron model. The black curve represents 
the non-linear gain function <j)(x) = 5 + 5 tanh(/3x). The dashed gray line is its tangent at 
the mean input value (denoted by the diagonal cross). The solid curve is the slope (<f)') 
averaged over the distribution of the fluctuating input (27). This distribution estimated 
from direct simulation is presented by black dots, the corresponding theoretical prediction 
of a normal distribution A/"(/U, <j) (27) is shown as the light gray curve. 



where 3i is the vector of incoming synaptic weights into neuron i and is a non-linear 
gain function. Assuming that the fluctuations of the total input J^n into the i-th neuron 
are sufficiently small to allow a linearization of the gain function <f>, we obtain the Taylor 
expansion 

F t (n,t) = F i (a)+cf)'(J i a)J i (xi(t-d)-a(t-d)) ) 

where 



*'(J<a) 



(26) 



is the slope of the gain function at the point of mean input. 

Up to this point the treatment of the system is identical to the work of Ginzburg & 
Sompolinsky (1994). Now we present an alternative approach for the linearization which 
takes into account the effect of fluctuations in the input. For sufficiently asynchronous 
network states, the fluctuations in the input J^n(t — d) to neuron i can be approximated 
by a Gaussian distribution A/"(/i, a). In the following we consider a homogeneous random 
network with fixed in-degree as described in "Population-averaged covarianccs" . As each 
neuron receives the same number K of excitatory and ~/K inhibitory synapses, the marginal 
statistics of the summed input to each neuron is identical. The mean input to a neuron 
then is /i = KJ{\ — jg)a, where a is the mean activity of a neuron in the network. If 
correlations are small, the variance of this input signal distribution can be approximated as 
the sum of the variances of the individual contributions from the incoming signals, resulting 
in a 2 = KJ 2 (1 + r yg 2 ) a(l — a), where we used the fact that the variance of a binary variable 
with mean a is a(l — a). This results from a direct calculation: since n € {0, 1}, n — n, so 
that the variance is (n 2 ) — (n) 2 = (n) — (n) 2 = a(l — a). Averaging the slope <f>' of the gain 
function over the distribution of the input variable results in the averaged slope 



A/"(/i,<T, x) 4> ' (x) dx 



(27) 



with A/*(/i, a, x) 



—00 

1 



27TCT 



exp 



2a 2 



The two alternative methods of linearization of <fi are illustrated in Figure 4. In the given 
example, the linearization procedure taking into account the fluctuations of the input signal 
results in a smaller effective slope (<//) than taking the slope </>'(a) at the mean activity a 
near its maximum. Averaging the slope (</>') over this distribution fits simulation results 
better than <p'(a) calculated at the mean of a, as shown in Figure 5. 

The finite slope of the non-linear gain function can be understood as resulting from the 
combination of a hard threshold with an intrinsic local source of noise. The inverse strength 
of this noise determines the slope parameter f3 (Ginzburg & Sompolinsky, 1994). In this 
sense, the network model contains two sources of noise, the explicit local noise, quantified 
by /3 and the fluctuating synaptic input interpreted as self-generated noise on the network 
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level, quantified by a. Even in the absence of local noise (/3 — > oo), the above mentioned 
linearization is applicable and yields a finite effective slope (</>') (27). In the latter case the 
resulting effective synaptic weight is independent of the original synapse strength (Grytskyy 
et al., 2013). 

We now extend the classical treatment of covariances in binary networks (Ginzburg 
& Sompolinsky, 1994) by synaptic conduction delays. In (25) Fj(n, £) must therefore be 
understood as a functional acting on the function n(i') for t' e [—00, t], so that also synaptic 
connections with time delay d can be realized. We define an effective weight vector to 
absorb the gain factor as w, = /3jJj, with either 0i = </>'(p) or /% = (</>') depending on the 
linearization procedure, and expand the right hand side of (24) to obtain 

N 

((nj(s)-aj(s))Fi(n,t)) = J~]wikCkj(t — d,s). 

Thus the cross-covariance fulfills the matrix delay differential equation 

r— c(t, s) +c(t,s) = wc(t-d,s). (28) 

at 

This differential equation is valid for t > s. For the stationary solution, the differential 
equation only depends on the relative timing u = t — s 

t—c(u)+c(u) = wc(u-d). (29) 

The same linearization applied to (23) results in the boundary condition for the solution of 
the previous equation 

2c(0) = wc(-d) + (wc(-d)) T (30) 

or, if we split c in its diagonal and its off-diagonal parts c a and c^ 

2c^(0) = wc^(-d) + (wc^(-d)) T + O (31) 

with O = wc (— d) + (wc a (— d)) T . 

In the following section we use this representation to demonstrate the equivalence of the 
covariance structure of binary networks to the solution for Ornstein-Uhlcnbcck processes 
with input noise. 

3.1. Equivalence of binary neurons and Ornstein-Uhlenbeck processes. In the 
following subsection we show that the same equations (29) and (31) for binary neurons also 
hold for the Ornstein-Uhlenbeck process (OUP) with input noise. In doing so here we also 
extend the existing framework of Ornstein-Uhlenbeck processes (Risken, 1996) to synaptic 
conduction delays d. A network of such processes is described by 

T-r(t)+r(t)=wr(i-rf)+x(£), (32) 

where x is a vector of pairwise uncorrelated white noise with (x(£)) x = and (xi(t)xj(t + 
t)) x = S i:j 5(t)p 2 . With the help of the Green's function G satisfying (rj- t + 1) G(t) = S(t), 
namely G(t) — - 9{t) e - */ T , we obtain the solution of equation (32) as 

r(t) = rG(t)r(0) + / G(t - t )(wr(* - d) + x(t )) dt . 
The equation for the fluctuations Sr(t) = r(t) — (r(t)) x around the expectation value 

8r(t)= / G(t~t )(wSr(t -d)+x(t))dt 
Jo 

coincides with the noisy rate model with input noise (4) with delay d and convolution kernel 
h = G. In the next step we investigate the covariance matrix Cij(t,s) — {5ri{t + s)Srj(t)) x 
to show for which choice of parameters the covariance matrices for the binary model and the 
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OUP with input noise coincide. To this end we derive the differential equation with respect 
to the time lag s for positive lags s > 

r*-c(t, s) = (r^Sr(t + s)Sr T (t)) x (33) 

as as 

= ((w5r(t + s-d)- Sr(t + s) + x(t + s))5r T (t)) x 
= wc(t, s — d) — c(t, s), 

where we used (x(< + s))5r(t)) x = 0, because the noise is realized independently for each 
time step and the system is causal. Eq. (33) is identical to the differential equation satisfied 
by the covariance matrix (28) for binary neurons (Ginzburg & Sompolinsky, 1994). To 
determine the initial condition of (33) we need to take the limit c(£, 0) = limg^+o c(i, s). 
This initial condition can be obtained as the stationary solution of the following differential 
equation 

rf t c(t, 0) = mn o ({rf t 6r(t + s)Sr T (t)) x + (Sr(t + s )r ±6r T (t)) x ) 

r T 

+ (Sr(t + s)(Sr T (t - d)w T - 5r T (t) + x T (t))) x ) 

= -2c(t, 0) + wc(i, -d) + c(t - d, d)w T + D. 

Here we used that (x(£ + s)<5r T (t)) vanishes due to independent noise realizations and causal- 
ity and 

D = lim (5r(t + s)x T (t)} x 

rt+s 

= lim / G{t + s-t')(w(5r(t - d)x(t)) x + (x(t')x(t)) x )dt 

s^+0, s<dj v v < v v / 

—0 causality = 15(t — t')p 2 

rt+s 



lim ({(wdr(t + s-d)- 6r(t + s) + x(t + s))5r T (t)} a 

S— ¥+0 



lim / G(t + s-t')l5(t-t')p 2 dt 
f0, s<dj 



s— >+0, s<d 

lim G(s)la 2 = -lp 2 . 

s->+0, s<d T 

In the stationary state, c only depends on the time lag s and is independent of the first time 
argument t, which, with the symmetry c(— d) T = c(d) yields the additional condition for 
the solution of (33) 

2c(0) = wc(-d) + (wc(-d)) T + D 
or, if c is split in diagonal and off-diagonal parts c a and c^, respectively, 

2c # (0) = wc^(-d) + {wc^(-d)) T + O 
2c a (0) = wcjti-d) + (wc # (-d)) T + D 

with O = wc a (- d) + (wc a (— d)) T . In the equation for the autocovariance c a the first two 
terms are contributions due to the cross covariance. In the state of asynchronous network 
activity with Cij ~ TV -1 for i ^ j these terms are typically negligible in comparison to the 
third term because ^ fc WikCki ~ wKN^ 1 = pw, which is typically smaller than 1 for small 
effective weights w < 1 and small connection probabilities p« 1. In this approximation 
with (33) the temporal shape of the autocovariance function is exponentially decaying with 
time constant r. With c o (0) « D/2 the approximate solution for the autocovariance is 

Co (t) = £exp(-M). (34) 

Z T 

The cross covariance then satisfies the initial condition 

2c # (0) = wc^(-ef) + (wc # (-d)) T + O 
O = wD/2 + (wD/2) T , 
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Figure 5: Binary model neuron corresponds to OUP model with input noise. Autoco- 
variance (A), crosscovarince (C), and autocovariance of population averaged activity (B,D) 
for binary neurons (dots) and rate model with input noise (crosses). C££,cgx and Cxx are 
shown in black, gray, and light gray. Corresponding theoretical predictions ((18) in B and 
D, (34) in A, their difference in C) are plotted as light gray curves throughout. Dashed 
curve in B represents the theoretical prediction using the linearization with the slope at the 
mean activity (26), the solid curve shows the results for the slope averaged over Gaussian 
distributed input fluctuations (27) . The spread of the simulation results for binary neurons 
in panel B result from different realizations of the random connectivity. 

which coincides with (31) for binary neurons if the diagonal matrix containing the zero 
time autocorrelations c a (0) for binary neurons is equal to D/2, i.e. if the amplitude of the 
input noise p 2 — 2ra(l — a) and the effective linear coupling satisfies w^ = /3jJj. Figure 5 
shows simulation results for population averaged covariance functions in binary networks 
and in networks of OUPs with input noise where the parameters of the OUP network are 
chosen according to the requirements derived above. The theoretical results (18) agree well 
with the direct simulations of both systems. For comparison, both methods of linearization, 
as explained above, are shown. The linearization procedure which takes into account the 
noise on the input side of the non-linear gain function results in a more accurate prediction. 
Moreover, the results derived here extend the classical theory (Ginzburg & Sompolinsky, 
1994) by considering synaptic conduction delays. Figure 7 shows the decomposition of the 
covariance structure for a non-zero delay d = 3 ms. For details of the implementation see 
"Implementation of binary neurons in a spiking simulator code" . 

4. Hawkes processes 

In the following section we show that to linear order the covariance functions in networks 
of Hawkes processes (Hawkes, 1971) are equivalent to those in the linear rate network with 
output noise. Hawkes processes generate spikes randomly with a time density given by 
r(i), where neuron i generates spikes at a rate r^(i), realized independently within each 
infinitesimal time step. Arriving spike trains s influence r according to 



r(i) =p+(h d *Js)(t), 



(35) 



with the connectivity matrix J and the kernel function hd including the delay. Here v is a 
constant base rate of spike emission assumed to be equal for each neuron. Here we employ 
the implementation of the Hawkes model in the NEST simulator (Gcwaltig & Diesmann, 
2007). The implementation is described in "Implementation of Hawkes neurons in a spiking 
simulator code". 
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Given neuron j spiked at time u < t, the probability of a spike in the interval [t, t + St) for 
neuron i is 1 if % = j, u = t (the neuron spikes synchronously with itself) and ri(t)St + o(St 2 ) 
otherwise. Considering the system in the stationary state with the time averaged activity 
f = (s(i)) we obtain a convolution equation for time lags r > for the covariance matrix 
with the entry c« (r) for the covariance between spike trains of neurons i and j 

c(r) = (s(i + r)s T (t)) - (s(i + r)>(s T (i)) (36) 

= ((,5(r)l + r(t + T))s T (t))-rf T 
= <r(t + r)(s T (t)-r T ))+D F 
= ((u + (h d * Js)(£ + r)(s T (t) - r T )) + D F 
= /i d *J(s(* + r)(s T (i)-r T ))+D F 
= (/i d *Jc)(r)+D F , 

with the diagonal matrix D F = <5(r)diag(r), which has been derived earlier (Hawkcs, 1971). 
If the rates of all neurons are equal, f.; = r, all entries in the diagonal matrix are the same, 
D r = (5(r)lf. In the subsequent section we demonstrate that the same convolution equation 
(36) holds for the linear rate with output noise. 

4.1. Convolution equation for linear noisy rate neurons. For the linear rate 
model with output noise we use equation (3) for time lags r > to obtain a convolution 
equation for the covariance matrix of the output signal vector y = r + x as 

c(r) = (y(t + r)(y T (t)-r T )) (37) 

= ((^*wy + x)(i + T)(y T (t)-f T )) 
= (h d * wc)(r) + (x(t + r)(r T (t) - r T )) + (x(t + r)x T (t)) 
= (ha *wc)(r) +D, 

where we utilized that due to causality the random noise signal generated at t + r has no 
influence on r(£), so the respective correlation vanishes. D is the covariance of the noise as 
in (11), Dij(r) = (xi(t)xj(t + r)) — 5ij5(t)p 2 . If p is chosen such that p 2 coincides with the 
averaged activity f in a network of Hawkes neurons and the connection matrix w is identical 
to J of the Hawkes network, the equations (36) and (37) are identical. Therefore the cross 
spectrum of both systems is given by (11). 

4.2. Non-linear self-consistent rate in rectifying Hawkes networks. The convo- 
lution equation (36) for the covariance matrix of Hawkes neurons is exact if no element of r is 
negative, which is particularly the case for a network of only excitatory neurons. Especially 
in networks including inhibitory couplings, the intensity r^ of neuron i may assume negative 
values. A neuron with r, < does not emit spikes, so the instantaneous rate is given by 
Aj = [r.;(£)] + = 0(ri(i))ri(t), with the Heaviside function 9. We now take into account this 
effective nonlinearity -the rectification of the Hawkes model neuron- in a similar manner 
as we already used to linearize binary neurons. If the network is in the regime of low spike 
rates, the fluctuations in the input of each neuron due to the Poissonian arrival of spikes are 
large compared to the fluctuations due to the time varying intensities r(t). Considering the 
same homogeneous network structure as described in "Population-averaged covariances" , 
the input statistics is identical for each cell i, so the mean activity Ao = (Aj) is the same 
for all neurons i. The superposition of the synaptic inputs to neuron i cause an instanta- 
neous intensity r, that follows approximately a Gaussian distribution M(p, a, r{) with mean 

p, = (r) = v + X Q KJ(1 - 57) and standard deviation a = yf{r 2 ) - (r) 2 = Jj^K(l + g 2 ~/). 

These expressions hold for the exponential kernel (5) due to Campbell's theorem (Papoulis 
& Pillai, 2002), because of the stochastic Poisson-like arrival of incoming spikes, where the 
standard deviation of the spike count is proportional to the square root of the intensity Aq. 
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The rate Ao is accessible by explicit integration over the Gaussian probability density as 
A = / Af(n,<r,r)r6(r)dr 



27TCT 



r=~ / eX P( T^2T- > rdr 

lira Jo 2cr^ 



2tt Jo ' 2cr 2 a 2 y^cr 7 2o- 



2^-^ 2cx 2 ' 2^ ~" v V2a' 

This equation needs to be solved self-consistently (numerically or graphically) to determine 
the rate in the network, as the right hand side depends on the rate Ao itself through \i and 
a. Rewritten as 

2 



1 1 c , l-i 

""" "71"a ; ' 



*V(r > 0) = - - -«*(—£-), (38) 



P/j,,a( r > 0) is the probability that the intensity of a neuron is above threshold and therefore 
contributes to the transmission of a small fluctuation in the input. A neuron for which r < 
acts as if it was absent. Hence we can treat the network with rectifying neurons completely 
analogous to the case of linear Hawkes processes, but multiply the synaptic weight J or —g J 
of each neuron with P^^ir > 0), i.e. the linearized connectivity matrix is 

w = P Mi<r (r>0)J. (39) 

Figure 6 shows the agreement of the covariance functions obtained from direct simulation 
of the network of Hawkes processes and the analytical solution (21) with average firing rate 
Ao determined by (38), setting the effective strength of the noise p 2 — A , and the linearized 
coupling as described above. The detailed procedure for choosing the parameters in the 
direct simulation is described together with the implementation of the Hawkes model in 
"Implementation of Hawkes neurons in a spiking simulator code" . 

5. Leaky integrate-and-fire neurons 

In this section we consider a network of leaky integrate-and-fire (LIF) model neurons with 
exponentially decaying postsynaptic currents and show its equivalence to the network of 
Ornstein-Uhlenbeck processes with output noise, valid in the asynchronous irregular regime. 
A spike sent by neuron j at time t arrives at the target neuron i after the synaptic delay 
d, elicits a synaptic current Ii that decays with time constant r s and causes a response 
in the membrane potential Vi proportional to the synaptic efficacy Jy. With the time 
constant T m of the membrane potential, the coupled set of differential equations governing 
the subthreshold dynamics of a single neuron i is (Fourcaud & Brunei, 2002) 

T m -^ - -Vi+Ii{t) 

at 

dl- N 

t s — j = -Ij + T m } j JijSjjt - d), (40) 

where the membrane resistance was absorbed into the definitions of J^- and 1^. If V^ reaches 
the threshold Vg at time point t\ the neuron emits an action potential and the membrane 
potential is reset to V r , where it is clamped for the refractory time r r . The spiking activity of 
neuron i is described by this sequence of action potentials, the spike train Si(t) = ^2 k S(t—t l k ). 
The dynamics of a single neuron is deterministic, but in network states of asynchronous, 
irregular activity and in the presence of external Poisson inputs to the network, the summed 
input to each cell can well be approximated as white noise (Brunei, 2000) with first moment 
Hi = T m J2j Jij r j an d second moment af = r m ^ ■ Jhrj, where Tj is the stationary firing 
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rate of neuron j. The stationary firing rate of neuron i is then given by (Fourcaud & Brunei, 
2002) 

r- 1 = T r + T m ^(F(y d )-F(y r )) (41) 

f(y) = e» 2 (l+erf(y)) F(y) = f f(y)dy 

with yg. r = 



« = V2|C(^)|, 



with Riemann's zeta function £. The response of the LIF neuron to the injection of an 
additional spike into afferent j determines the impulse response Wijh(t) of the system. The 
time integral Wij — Wij L h(t) dt is the DC-susceptibility, which can formally be written as 
the derivative of the stationary firing rate by the rate of the afferent Tj , which, evaluated by 
help of (41), yields (Helias et al., 2013, Results and App. A) 



Wij = ^ = aJij + fiJl (42) 



with a = \fTr{T m r l ) 2 —{f{y e )-f{y r )) 
and j3 = Atw,) 2 - -j [ f{ye) — - - f{y r ) — '- 



2aj 

In the strongly fluctuation-driven regime, the temporal behavior of the kernel h is domi- 
nated by a single exponential decay, whose time constant can be determined empirically. 
In a homogeneous random network the firing rates of all neurons arc identical n = f and 
follow from the numerical solution of the self-consistency equation (41). Approximating the 
autocovariance function of a single spike train by a J-peak scaled by the rate fS(t), one 
obtains for the covariance function c between pairs of spike trains the same convolution 
equation (36) as for Hawkes neurons (Helias ct al., 2013, cf. eq. 5). As shown in "Con- 
volution equation for linear noisy rate neurons" this convolution equation coincides with 
that of a linear rate model with output noise (37), where the diagonal elements of D are 
chosen to agree to the average spike rate p 2 — f. The good agreement of the analytical cross 
covariance functions (21) for the OUP with output noise and direct simulation results for 
LIF are shown in Figure 6. 

6. Discussion 

In this work we present a unified theoretical view on pairwise correlations in recurrent 
networks. We consider binary neuron models, leaky integrate-and-fire models, and linear 
point process models. These models containing a non-linearity (spiking threshold in spik- 
ing models, non-linear sigmoidal gain function in binary neurons, strictly positive rates in 
Hawkes processes) are linearized, taking into account the distribution of the fluctuating in- 
put, thereby extending the classical procedure which neglects these fluctuations. For the 
example of binary neurons, we show that capturing the influence of fluctuations in the 
effective system parameters explains simulation results better and moreover precisely deter- 
mines the onset and parameters of oscillations. A similar averaging approach allows us to 
extend the theory of correlations in networks of Hawkes neurons to inhibition. The novel 
approach enables us to demonstrate the equivalence of each of these models after linear 
approximation to a linear model with fluctuating continuous variables. For the commonly 
appearing exponentially decaying response kernel function, these rate models coincide with 
the Ornstein-Uhlenbeck process (OUP, Uhlcnbcck & Ornstein, 1930; Riskcn, 1996). We 
find that the considered models form two groups, which, in linear approximation merely 
differ by a matrix valued factor scaling the noise and in the choice of variables interpreted as 
neural activity. The difference between these two groups corresponds to the location of the 
noise: spiking models -leaky integrate-and-fire models and Hawkes models- belong to the 
class with noise on the output side, added to the activity of each neuron. The non-spiking 
binary neuron model corresponds to an OUP where the noise is added on the input side of 
each neuron. The closed solution for the correlation structure of OUP holds for both classes. 
We identify different contributions to correlations in recurrent networks: the solution for 
output noise is split into three terms corresponding to the <5-peak in the autocovariance, the 
covariance caused by shared input, and the direct synaptic influence of stochastic fluctuations 
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Figure 6: Covariance structure in spiking networks corresponds to OUP with output noise. 
A Autocovariance obtained by direct simulation of the LIF (black), Hawkes (gray), and OUP 
(light gray) models for excitatory (dots) and inhibitory neurons (crosses) . B Covariance csx 
averaged over disjoint pairs of neurons for LIF (black dots), Hawkes (gray dots), and OUP 
with output noise (empty circles). C Covariance averaged over disjoint pairs of neurons of 
the same type. D Autocovariance of the population averaged activity. Averages in C,D 
over excitatory neurons as black dots, over inhibitory neurons as gray dots. Corresponding 
theoretical predictions (21) are plotted as light gray curves in all panels except A. Light gray 
diagonal crosses in A and D denote theoretical peak positions determined by the firing rate 
f as fAt (where At = 0.1 ms is the time resolution of the histogram). 
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Figure 7: Different echo terms for spiking and non-spiking neurons. Binary non-spiking 
neurons shown in A,C and LIF in B,D. A,B Echo terms by direct influence of the neuron's 
output on the network in dependence of neuron types (in A, B C££,C£i, and c%x are plotted 
as black, gray dots and circles). C,D Contributions to the covariance evoked by correlated 
and common input (black dots) measured with help of auxiliary model neurons which do 
not provide feedback to the network. Corresponding theoretical predictions (16) are plotted 
as light gray curves throughout. 

of one neuron on another - the latter echo terms are equal to propagators acting with delays 
(Helias et al., 2013). A similar splitting into echo and correlated input terms for the case of 
input noise is shown in Figure 7. For increasing network size N — > oo, keeping the connection 
probability p fixed, so that K — pN, and with rescaled synaptic amplitudes J ~ 1/vN (van 
Vreeswijk & Sompolinsky, 1996; Rcnart ct al., 2010) the echo terms vanish fastest. Formally 
this can be seen from (18): the multiplicative factor of the common covariance term ip^ does 
not change with N while the other coefficients decrease. So ultimately all four entries of 
the matrix c have the same time dependence determined by the common covariance term 
if4,. In particular the covariance between excitation and inhibition csx becomes symmetric 
in this limit. This finally provides a quantitative explanation of the observation made in 
(Rcnart et al., 2010) that the time-lag between excitation and inhibition vanishes in the 
limit of infinitely large networks. For a different synaptic rescaling J ~ N~ l while keeping 
p 2 constant by appropriate additional input to each neuron (see Helias ct al., 2013, applied 
to the LIF model), all multiplicative factors decrease ~ N^ 1 and so does the amplitude of 
all covariances. Hence the asymmetry of csx does not vanish in this limit. The same results 
hold for the case of output noise where the term with tpi describes the common input part 
of the covariance. In this case and for finite network size, Cxs coincides with ess and cgx 
with cxx for t > 0, having a discontinuous jump at the time of the synaptic delay t = d. For 
time lags smaller than the delay all four covariances coincide. This is due to causality, as the 
second neuron cannot feel the influence of a fluctuation that happened in the first neuron 
less than one synaptic delay before. The covariance functions for systems corresponding to 
an OUP with input noise contain neither discontinuities nor sharp peaks at t = d, but Csx 
and cxs have maxima and minima near this location. This observation can be interpreted 
as a result of the stochastic nature of the binary model where changes in the input influence 
the state of the neuron only with a certain probability. So, the entries of c in this case take 
different values for |t| < d but show the tendency to approach each other with increasing 
\t\ 3> d. This tendency increases with network size. Our analytical solutions (18) for input 
noise and (21) for output noise hence explain the model-class dependent differences in the 
shape of covariance functions. 

Our approach enables us to map results obtained for one neuron model to another, 
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in particular we extend the theory of all considered models to capture synaptic conduction 
delays, and devise a simpler way to obtain solutions for systems considered earlier (Ginzburg 
& Sompolinsky, 1994). Our derivation of covariances in spiking networks does not rely on 
the advanced Wiener-Hopf method (Hazcwinkel, 2002), as earlier derivations (Hawkes, 1971; 
Helias ct al., 2013) do, but only employs elementary methods. Our results are applicable 
for general connectivity matrices, and for the purpose of comparison with simulations we 
explicitly derive population averaged results. The averages of the dynamics of the linear 
rate model equations are exact for random network architectures with fixed out-degree, and 
approximate for fixed in-degree. Still, for non-linear models the linearization for fixed in- 
degree networks are simpler, because the homogeneous input statistics results in an identical 
linear response kernel for all cells. Finally we show that the oscillatory properties of networks 
of integrate-and-fire models (Brunei, 2000; Helias et al., 2013) are model-invariant features 
of all of the studied dynamics, given inhibition acts with a synaptic delay. We relate the 
collective oscillations to the pole structure of the cross spectrum, which also determines the 
power spectra of population signals such as EEG, ECoG, and the LFP. 

The presented results provide a further step to understand the shape and to unify the 
description of correlations in recurrent networks. We hope that our analytical results will 
be useful to constrain the inverse problem of determining the synaptic connectivity given 
the correlation structure of neurophysiological activity measurements. Moreover the ex- 
plicit expressions for covariance functions in the time domain are a necessary prerequisite 
to understand the evolution of synaptic amplitudes in systems with spike-timing dependent 
plasticity and extend the existing methods (Burkitt et al., 2007; Gilson et al., 2009, 2010) 
to networks including inhibitory neurons and synaptic conduction delays. 
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8. Appendix 

8.1. Explicit calculation of the expression for the population averaged cross 
covariance in time domain. We obtain the population averaged cross spectrum for the 
Ornstein-Uhlenbeck process with input noise by inserting the averaged connectivity matrix 
w = M (14) into (8). The two eigenvalues of M are and L = Kw(l — 73). Taking these 
into account, we first rewrite the term 



(H^lu)- 1 -M)- 1 
= det(.ff d (u;)- 1 -M)- 1 



Hd(uj) 1 + Kwjg —Kwjg 

Kw Hdito)- 1 - Kw 



where we introduced f(to) — (iJ^cj) -1 — L)^ 1 . The corresponding transposed and conjugate 
complex term follows analogously. Hence we obtain the expression for the cross spectrum 
(17). The residue of J(uj) at cj = Zk{L) is 

Res(/,cj = z k (L)) = lim — - 



Wl->U) / l( Wl ) 

1'Hopitai .. 1 (d(e iud (l + iu>T)) 

lim 



ui-*u (/ 1 )'(wi) \ dio 

= (ide iujd {l + Iujt) + ire lud y l = (idL + ire iujd y l , 

where in the last step we used the condition for a pole Hd(zk)^ 1 — e lZkd (l + izut) = L (see 
"Spectrum of the dynamics"). The residue of Hd{w) at z(0) = - is —-e d l T . Using the residue 
theorem, we need to sum over all poles within the integration contour {zk{L)\k G N} U - 
to get the expression for c(t) = ^ f^ C(uj)e wt duj = iT, ze {z k (L)\keN}u± Res(C(z), z)e tzt 
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for t ^ 0. Sorting (17) to obtain four matrix prefactors and remainders with different 
frequency dependence, $i(w) = /(w)/(-w), $ 2 (w) = f{ui)f(-uj)H d (uj), $ 3 (w) = $ 2 (-w), 
and $4(0;) = f(u))f(—w)Hd(<u)Hd(—u), we get (18). C(w) for output noise (20) is obtained 
by multiplying the expression for C(w) for input noise with Hj(u))Hj(—u) = (1 + uj 2 t 2 ). 
In order to perform the back Fourier transformation one first needs to rewrite the cross 
spectrum in order to isolate the frequency independent term and the two terms that vanish 
for either t < d or t > d, as described in "Fourier back transformation" , 

C(«) =/(w)(l + Kw(™ ~2l ) H d ( W ))MDM T /(- W )(l + Kw f 2 9 ig \ ) H d {-u)) 

+/(«)(1 +Kw(™ ~2l ) H d (Lo))MT> 

+BM T f(-oj)(l + Kw( 2 9 ig \ ) H d (-uj)) + D 

=/(u)MDM T /(-u) + /(cj)MD + DM T /(-w) + D, 

/ ^yo — 'yq \ 
where in the last step we used ' I M = 0, because M is symmetric, obtaining 

(21). For each of the first three terms in the last expression the right integration contour 
needs to be chosen as described in "Fourier back transformation" on the example of the 
general expression (16). 

8.2. Implementation of noisy rate models. The dynamics is propagated in time 
steps of duration At (note that in other works we use ft, as a symbol for the computation 
step size, which here is used as the symbol for the kernel). The product of the connectivity 
matrix with the vector of output variables at the end of the previous step z — 1 is the vector 
I(ii) of inputs at the current step i. The intrinsic time scale of the system is determined by 
the time constant r. For sufficiently small time steps At<Cr these inputs can be assumed 
to be time independent within one step. So we can use (3) or (4) and analytically convolve 
the kernel function h assuming the input to be constant over the time interval At. This 
corresponds to the method of exponential integration (Rotter & Diesmann, 1999, see App. 
C.6) requiring only local knowledge of the connectivity matrix w. Note that this procedure 
becomes exact for At — > and for finite At is an approximation. The propagation of the 
initial value rj(tj-i) until the end of the time interval takes the form r.,(tj_i) e~ Af / r because 
h(ti) — h(ti-i) e~ At / r , so we obtain the expression rj{ti) at the end of the step as 

r 3 {U) = e~ At/r r^ti-i) + (1 - e" At/T ) /,-(*,), (43) 

where Ij denotes the input to the neuron j. For output noise the output variable of neuron 
j is yj = Tj +Xj, with the locally generated additive noise Xj and hence the input is Ij(t%) = 
(wy(tj))j. In the case of input noise the output variable is Tj and the additional noise is 
added to the input variable, Ij(U) = (wr(ij))j + Xj(ti). In both cases Xj is implemented 
as a binary noise: in each time step, Xj is independently and randomly chosen to be 1 or 
— 1 with probability 0.5 multiplied with p/yAt to satisfy (2) for discretized time. Here the 
<5-function is replaced by a "rectangle" function that is constant on the interval of length At, 
vanishes elsewhere, and has unit integral. The factor At -1 in the expression for x 1 ensures 
the integral to be unity. So far, the implementation assumes the synaptic delay to be zero. 
To implement a non-zero synaptic delay d, each object representing a neuron contains an 
array b of length l d = dj At acting as a ring buffer. The input Ij(U) used to calculate the 
output rate at step i according to (43) is then taken from position i mod l d of this array 
and after that replaced by the input presently received from the network, so that the new 
input will be used only after one delay has passed. This sequence of buffer handling can be 
represented as 




for input noise 
for output noise 

The model is implemented in Python version 2.7 (Python Software Foundation, 2008) using 
numpy 1.6.1 (Ascher et al., 2001) and scipy 0.9.0 (Jones et al., 2001). 
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8.3. Implementation of binary neurons in a spiking simulator code. The binary 
neuron model is implemented in the NEST simulator, version 2.2.1 (Gcwaltig & Diesmann, 
2007), which allows distributed simulation on parallel machines and handles synaptic delays 
in the established framework for spiking neurons (Morrison et al., 2005). The name of 
the model is "ginzburg_neuron" . In NEST information is transmitted in form of point 
events, which in case of binary neurons are sent if the state of the neuron changes: one 
spike is sent for a down-transition and two spikes at the same time for an up-transition, 
so the multiplicity reflects the type of event. The logic to decode the original transitions 
is implemented in the function handle shown in Alg. 2. If a single spike is received, the 
synaptic weight w is subtracted from the input buffer at the position determined by the 
time point of the transition and the synaptic delay. In distributed simulations a single spike 
with multiplicity 2 sent to another machine is handled on the receiving side as two separate 
events with multiplicity 1 each. In order to decode this case on the receiving machine we 
memorize the time (ti ast ) and origin (global id gid last of the sending neuron) of the last 
arrived spike. If both coincide to the spike under consideration, the sending neuron has 
performed an up transition — > 1. We hence add twice the synaptic weight 2w to the 
input buffer of the target neuron, one that reflects the real change of the system state and 
another that compensates the subtraction of w after reception of the first spike of a pair. 
The algorithm relies on the fact that within NEST two spikes that are generated by one 
neuron at the same time point are delivered sequentially to the target neurons. This is 
assured, because neurons are updated one by one: The update propagates each neuron by a 
time step equal to the minimal delay d m i n in the network. All spikes generated within one 
update step are written sequentially into the communication buffers, and finally the buffers 
are shipped to the other processors (Morrison et al., 2005). Hence a pair of spikes generated 
by one neuron within a single update step will be delivered consecutively and will not be 
interspersed by spikes from other neurons with the same time stamp. 

The model exhibits stochastic transitions (at random points in time) between two states. 
The transitions are governed by probabilities cj>(h). Using asynchronous update (Rumelhart 
et al., 1986), in each infinitesimal interval [t,t + St) each neuron in the network has the 
probability -St to be chosen for update (Hopfield, 1982). A mathematically equivalent 
formulation draws the time points of update independently for all neurons. For a particular 
neuron, the sequence of update points has exponentially distributed intervals with mean 
duration r, i.e. it forms a Poisson process with rate r _1 . We employ the latter formulation 
to incorporate binary neuron models in the globally time-driven spiking simulator NEST 
(Gewaltig & Diesmann, 2007) and constrain the points of transition to a discrete time grid 
At = 0.1 ms covering the interval d m in _ At. This neuron state update is implemented by 
the algorithm shown in Alg. 1. Note that the field h is updated in steps of At while the 
activity state is updated only when the current time exceeds the next potential transition 
point. As the last step of the activity update we draw an exponentially distributed time 
interval to determine the new potential transition time. The potential transition time is 
represented with a higher resolution (on the order of microseconds) than At to avoid a 
systematic bias of the mean inter-update-interval. This update scheme is identical to the 
one used in (Hopfield, 1982). Note that the implementation is different from the classical 
asynchronous update scheme (van Vreeswijk & Sompolinsky, 1998), where in each discrete 
time step At exactly one neuron is picked at random. The mean inter-update-interval (time 
constant r in Alg. 1) in the latter scheme is determined by r = AtN, with N the number of 
neurons in the network. For small time steps both schemes converge so that update times 
follow a Poisson process. 

At each update time point the neuron state becomes 1 with the probability given by 
the function <f> applied to the input at that time according to (25) and with probability 
1 — 0. The input is a function of the whole system state and is constant between spikes 
which indicate state changes. Each neuron therefore maintains a state variable h at each 
point in time holding the summed input and being updated by adding and subtracting the 
input read from the ring buffer b at the point readpos(t) corresponding to the current time 
(see Morrison et al., 2005, for the implementation of the ring buffer, i.p. Fig 6). The ring 
buffer enables us to implement synaptic delays. For technical reasons this implementation 
requires a minimal delay of a single simulation time step (Morrison & Diesmann, 2008). The 
gain function (j) applied to the input h has the form 

4>(h) = Cl h + c 2 - (1 + tanh(c 3 (/i - 0))), (44) 
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where throughout this manuscript we used ci = 0, eg = 1, and C3 = /?, as defined in 
"Parameters of simulations" . 

Algorithm 1 Update function of a binary neuron embedded in the spiking network simulator 
NEST. The function readpos(t) returns a position in the ring buffer 6 corresponding to the 
current time point. 

1 y 4— // initially neuron is inactive 

2 tnext < rlog(rand()) // next time point of update 

3 

4 for each time step t: 

5 

6 h <— h + 6[readpos(i)] 

7 

8 i f t > inoxt : 

9 // up— state with probability given by 

10 // gain function <f> depending on input h(t) 

11 

12 if 4>{h) >rand(): 

13 J/ncw <~ 1 

14 else : 

15 y ncw <- 
16 

17 i f y„ow ^ V : 

18 // down transition : send single spike 

19 // up transition: send two spikes 

20 send (y now + 1 spikes) 
21 

22 y <r- y nC w 
23 

24 // add an exponentially distributed time interval 

25 i„oxt <- ^ncxt - t log(rand()) 



8.4. Implementation of Hawkes neurons in a spiking simulator code. Hawkes 
neurons (Hawkes, 1971) were introduced in the NEST simulator in version 2.2.0 (Gewaltig 
& Diesmann, 2007). The name of the model is "pp_psc_delta". In the following we describe 
the implemented neuron model in general and mention the particular choices of parameter 
and correspondences to the theory presented in "Hawkes processes" . The dynamics of the 
quasi-membrane potential u is integrated exactly within a time step At of the simulation 
(Rotter & Diesmann, 1999), expressing the voltage u(tj) at the end of time step i by the 
membrane potential at the end of the previous time step u(tj_i) as 

u{U) = e- At / r u(*i_i) + (1 - e' At / r ) R m I e + b(U), (45) 

where I e is a time-step wise constant input current (equal to in all simulations presented 
in this article) and R m = T m /C m is the membrane resistance. The buffer b(ti) contains the 
summed contributions of incoming spikes, multiplied by their respective synaptic weight, 
which have arrived at the neuron within the interval (t i _ 1 ,ij]. b is implemented as a ring- 
buffer in order to handle the synaptic delay, logically similar as in "Implementation of 
noisy rate models", described in detail in Morrison et al. (2005). The instantaneous spike 
emission rate is A = [c\U + C2e C3 "] + , where we use C3 = in all simulations presented 
here. The quantities in the theory "Hawkes processes", in particular in (35), are related to 
the parameters of the simulated model in the following way The quantity r relates to the 
membrane potential tiasr = C1U + C2 and the background rate v agrees to C2 = v. Hence the 
synaptic weight J^- corresponds to the synaptic weight in the simulation multiplied by c\. 
For the correspondence of the Hawkes model to the OUP with output noise of variance p 2 
we use (38) to adjust the background rate v in order to obtain the desired rate Ao = p 2 and 
we choose the synaptic weight J of the Hawkes model so that the linear coupling strength 
w of the OUP agrees to the effective linear weight given by (39). These two constraints 
can be fulfilled simultaneously by solving (38) and (39) by numerical iteration. The spike 
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Algorithm 2 Input spike handler of a binary neuron embedded in the spiking network 
simulator NEST. The simulation kernel calls the handle function for each spike event to 
be delivered to the neuron. A spike event is characterized by the time point of occurrence 
^spikei the synaptic delay d after which the event should reach the target, the global id 
gid identifying the sending neuron, and the multiplicity m > 1, indicating the reception of 
multiple spike events. The function pos(£ sp ikc, d, t) returns the position in the ring buffer b 
to which the spike is added so that it will be read at time t + d by the update function of 
the neuron, see Alg. 1. 

1 handle (t spikc ,d, gid, m) : 

2 

3 if m=l: 

4 

5 // multiplicity = 1, either a single 1 — » event 

6 // or the first or second of a pair of — >• 1 events 
7 

8 i f gid = gid lastspikG and i spiko = ii as ts P ike : 
9 

10 // received twice the same event, so transition — >• 1 

11 // add 2w to compensate for subtraction after reception 

12 // of first event 

13 6[pos(t sp iko, d, t)] <- 6[pos(i sp ike, d, t)] + 2w 
14 

15 else : 

16 // count this event negatively , 

17 // assuming it comes as single event 

18 // transition 1 — ^ 

19 6[pos(t sp ikc, d, t)] i- 6[pos(i spi ke, d, t)] - w 

20 else : 
21 

22 // multiplicity != 1 

23 

24 i f m = 2 : 

25 

26 // count this event positively , transition — >• 1 

27 6[pos(i spi kc, d, t)} «- %os(£ sp ike, d, t)} + w 

28 gidiastspikc <~ g id 

^y Mastspikc ^ ^spikc 



emission of the model is realized either with or without dead time. In this article we only 
used the latter. In the presence of a dead time, which is constrained to be larger than 
the simulation time step, at most one spike can be generated within a time step. A spike 
is hence emitted with the probability p>i = 1 — e AAt , where e AAt is the probability of the 
complementary event (emitting spikes) , implemented by comparing a uniformly distributed 
random number to p>\. The refractory period is handled as described in Morrison et al. 
(2005). Without refractoriness, the number of emitted spikes is drawn from a Poisson 
distribution with parameter A At, implemented in the GNU Scientific Library (Galassi et al., 
2006). Reproducibility of the random sequences for different numbers of processes and 
threads is ensured by the concept of random number generators assigned to virtual processes, 
as described in (Plesser et al., 2007). 

8.5. Parameters of simulations. For all simulations we used 7 = 0.25 corresponding 
to the biologically realistic fraction of inhibitory neurons, a connectivity probability p = 0.1, 
and a simulation time step of At = 0.1 ms. For binary neurons we measured the correlation 
functions with a resolution of 1 ms, for all other models the resolution is 0.1 ms. The 
covariance is obtained for a time window of ±100 ms. 

The parameters for simulations of the LIF model presented in Figure 6 and Figure 7 are 
J = 0.1 mV, t = 20ms, t s = 2ms, r r = 2ms, V e = 15mV, V r = 0, g = 6, d = 3ms, N = 8000. 
The number of neurons in the corresponding networks of other models is the same. Cross 
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covariances are measured between the summed spike trains of two disjoint populations of 
N TCC = 1000 neurons each. The single neuron autocovariances a a are averaged over a 
subpopulation of 100 neurons. The autocovariances of the population averaged activity 
jj-a a +C aa for population a € {S,X} (shown in Figure 6) are constructed from the estimated 
single neuron population averaged autocovariances a a and cross covariances C aa . This 
enables us to estimate a a and C aa from the activity of a small subpopulation and still assigns 
the correct relative weights to both contributions. The corresponding effective parameters 
describing the system dynamics are p, — 15 mV, a = 10 mV, r — 23.6 Hz (see (40) and the 
following text for details). 

The parameters of the Hawkes model and of the noisy rate model with output noise 
yielding quantitatively agreeing covariance functions are: 

• For simulations of the noisy rate model with output noise presented in Figure 6 and 
Figure 2 the parameters are w = 0.0043, g « 5.93, r = 4.07ms, p 2 — 23.6 Hz, d = 3ms 
(see (3), (4)). In Figure 2 also results for d = 1 ms and for input noise are shown. 
Signals are measured from 7V rec = 500 neurons in each population to obtain cgx, cxs 
and from the whole population to determine egg and cxx- The cross covariances C$e 
and Cxx are estimated from two disjoint subpopulations each comprising half of the 
neurons of the respective population. 

• For the network of Hawkes neurons presented in Figure 6 we used A « 22.54 Hz (see 
(38)), J = 0.0055 mV, d — 3 ms, and the same g and r as for the noisy rate model. 
We measured the cross covariances in the same way as for the LIF model, but using 
the spike trains from sub-populations of N TSC — 2000 neurons. The autocovariances of 
the population averaged activity were estimated from the whole populations. 

The network of binary neurons shown in Figure 7 uses 6 = —3.89 mV, /3 = 0.5 mV , 
J = 0.02 mV, d — 3 ms (see (25), (44)), and the same g and r as the noisy rate model. 
Covariances are measured using the signals from all neurons. 

The simulation results for the network of binary neurons presented in Figure 5 uses 
9 = -2.5 mV, t = 10 ms, /3 = 0.5 mV -1 , g = 6, J « 0.0447 mV, N = 2000 and the smallest 
possible value of synaptic delay is d — 0.1 ms equal to time resolution (the same set of 
parameters only with modified /3 = 1 mV~ was used to create Figure 4). The cross covari- 
ances Css an( i Cxi are estimated from two disjoint subpopulations each comprising half of 
the neurons of the respective population, csx is measured between two such subpopulations. 
For C££ and cxx we used the full populations. 

The parameters required for a quantitative agreement with the rate model with input 
noise are w ~ 0.011, p « 2.23 y'ms. We used the same parameters in Figure 3, where 
additionally results for w — 0.018 are shown. The population sizes are the same as for the 
binary network. The covariances are estimated in the same way as for the rate model with 
output noise. Note that the definition of noisy rate models has no limitation for units of p 2 . 
These can be arbitrary and are chosen differently as required by the correspondence with 
either spiking or binary neurons. 
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